Single-cell RNA sequencing data processing

In this example, we are using a CITEseq dataset of mouse glioblastoma model.This dataset contains CD45+ sorted cells from a GL261 glioblastoma tumor of wild type and CCR2 knockout mice (WT4 and KO4).You can download the expression matrix here.

suppressPackageStartupMessages({
  library(Seurat)
  library(scater)
  library(dplyr)
  library(VennDiagram)
  library(mvoutlier)
  library(harmony)
  library(clustree)
  library(cowplot)
  library(scDblFinder)
  library(ggplot2)
  library(ggrepel)
})

Set the working directory

knitr::opts_knit$set(root.dir = "path/to/directory")

1) LOADING THE INPUT DATA

We use the Read10X function to read in the expresion matrix from the Cell Ranger output, which returns a list of two matrices: “Gene Expression” and “Antibody Capture”. We load the “Gene Expression” matrix into a Seurat object.

seur<-CreateSeuratObject(counts = Read10X(paste0("data/Citeseq_mouse_GBM/filtered_feature_bc_matrix"))[["Gene Expression"]] )

10X data contains more than one type and is being returned as a list containing matrices of each type.

### For single-cell RNA sequencing experiment wihout feature barcoding, use:
### seur<-CreateSeuratObject(counts = Read10X(paste0("path/to/filtered_feature_bc_matrix")) )

Identify the different samples based on the cell barcode names

head(colnames(seur))

[1] "AAACCCACAAGCGCTC-1" "AAACCCACAAGTCCCG-1" "AAACCCACACTGGCCA-1" "AAACCCACAGAGGACT-1" "AAACCCACAGCTGGTC-1" "AAACCCACAGTAGAAT-1"

seur$sample<-sapply(strsplit(colnames(seur), "-"), "[[",2)
table(seur$sample)

    1     2 
15850 12434 

seur$sample<-plyr::mapvalues(seur$sample, from=1:2, to=c("KO1","WT1"))

Alternatively, if starting from the individual expression matrices for each sample,we can load them as a list of seurat objects and then merge them using the Seurat function “merge”:

# sample_folder_names=c("KO1","WT1")
### sample_folder_names - a vector with the names of the folders that contain the Cellranger outputs for each sample
# seur=list()
# for ( i in c(sample_folder_names)) {  
#   seur[[i]]<-CreateSeuratObject(counts = Read10X(paste0(i,"/filtered_feature_bc_matrix"))[["Gene Expression"]] )
#   seur[[i]]$sample=i
# }
# seur<-merge(seur[[1]],seur[2:length(seur)])

2) QC: CELLS

2.1) Find outliers for total UMI counts, number of genes and percent mitohondrial genes per cell.

Outliers are determined based on the median absolute deviation (MAD). The goal is to remove the long tail before/after the peak in the sitribution of the QC metric. It is possible to modify the nmads parameter (minimum number of MADs away from median, required for a value to be called an outlier), or to set the threshold manually (e.g. remove all cells with percent mitochondrial genes above 40).

Calculate % mitochondrial genes per cell

seur[["perc.mito"]] <- PercentageFeatureSet(seur, pattern = "^mt-")
summary(seur$perc.mito)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   3.179   4.701   7.913   6.769  96.575 

- UMI counts per cell

outliers=c()
for ( i in unique(seur$sample)){
  outliers=c(outliers, 
             scater::isOutlier(seur$nCount_RNA[seur$sample==i], nmads=3, type="lower", log=TRUE)
             )
}
seur$nUMI.outlier.low <- outliers[colnames(seur)]
cat("Outliers:",sum(seur$nUMI.outlier.low))

Outliers: 1

Create histograms

for ( i in unique(seur$sample)){
  hist(seur$nCount_RNA[seur$sample==i],
        breaks = 100,xlab="nCount_RNA",
        main=paste0("Total UMI counts per cell: ",i))
  if(sum(seur$sample==i & seur$nUMI.outlier.low)!=0)
    abline(v = max(seur$nCount_RNA[seur$sample==i & seur$nUMI.outlier.low]), col = "red")
}

Create violin plots

for ( i in unique(seur$sample)){
 print(ggplot(as.data.frame(seur[[]])[seur$sample==i,], aes(1, nCount_RNA)) + 
    geom_violin(fill="gray80") +theme_classic()+ theme(axis.title.x = element_blank())+
    geom_jitter(height = 0, width = 0.3, aes(col=nUMI.outlier.low)) +
    scale_color_manual(values=c("#00BFC4", "#F8766D"))+ggtitle(paste0("Total UMI counts per cell: ",i)))
}

- Number of genes per cell

outliers=c()
for ( i in unique(seur$sample)){
  outliers=c(outliers, 
             scater::isOutlier(seur$nFeature_RNA[seur$sample==i], nmads=3, type="lower", log=TRUE)
             )
}
seur$nGene.outlier.low <- outliers[colnames(seur)]
cat("Outliers:",sum(seur$nGene.outlier.low))

Outliers: 1009

Create histograms

for ( i in unique(seur$sample)){
  hist(seur$nFeature_RNA[seur$sample==i],
        breaks = 100,xlab="nCount_RNA",
        main=paste0("Number of genes per cell: ",i))
  if(sum(seur$sample==i & seur$nGene.outlier.low)!=0)
    abline(v = max(seur$nFeature_RNA[seur$sample==i & seur$nGene.outlier.low]), col = "red")
}

Create violin plots

for ( i in unique(seur$sample)){
 print(ggplot(as.data.frame(seur[[]])[seur$sample==i,], aes(1, nFeature_RNA)) + 
    geom_violin(fill="gray80") +theme_classic()+ theme(axis.title.x = element_blank())+
    geom_jitter(height = 0, width = 0.3, aes(col=nGene.outlier.low)) +
    scale_color_manual(values=c("#00BFC4", "#F8766D"))+ggtitle(paste0("Number of genes per cell: ",i)))
}

- Proportion of mitohondrial genes per cell

outliers=c()
for ( i in unique(seur$sample)){
  outliers=c(outliers, 
             scater::isOutlier(seur$perc.mito[seur$sample==i], nmads=3, type="higher", log=TRUE)
             )
}
seur$mito.outlier.high <- outliers[colnames(seur)]
cat("Outliers:",sum(seur$mito.outlier.high))

Outliers: 1440

Create histograms

for ( i in unique(seur$sample)){
  hist(seur$perc.mito[seur$sample==i],
        breaks = 100,xlab="perc.mito",
        main=paste0("% mito genes per cell: ",i))
  if(sum(seur$sample==i & seur$mito.outlier.high)!=0)
    abline(v = min(seur$perc.mito[seur$sample==i & seur$mito.outlier.high]), col = "red")
}

Create violin plots

for ( i in unique(seur$sample)){
 print(ggplot(as.data.frame(seur[[]])[seur$sample==i,], aes(1, perc.mito)) + 
    geom_violin(fill="gray80") +theme_classic()+ theme(axis.title.x = element_blank())+
    geom_jitter(height = 0, width = 0.3, aes(col=mito.outlier.high)) +
    scale_color_manual(values=c("#00BFC4", "#F8766D"))+ggtitle(paste0("% mito genes per cell: ",i)))
}

Overlap of cells, outliers for UMI counts, number of genes and % mitochodrial genes per cell

A large nuber of cells, outliers specifically for one QC metric might be a concern and require fuurhter investigation

for ( i in unique(seur$sample)){
v <-venn.diagram(
  list (nUMI=colnames(seur)[seur$nUMI.outlier.low & seur$sample==i],
        nGene=colnames(seur)[seur$nGene.outlier.low & seur$sample==i],
    perc.mito=colnames(seur)[seur$mito.outlier.high & seur$sample==i]),
  filename=NULL,main=i,
  alpha = c( 0.5,0.5,0.5),
  fill = c("green","orange","blue")
)
grid.newpage()
grid.draw(v)
rm(v)
}

How many cells will be removed during the filtering

cells.to.keep= ! (seur$mito.outlier.high | seur$nGene.outlier.low | seur$nUMI.outlier.low)
print(paste(ncol(seur)- 
            ncol(seur[,cells.to.keep]),
     "cells to be removed"))

[1] "1524 cells to be removed"

Violin plots after filtering

for ( i in unique(seur$sample)){
  print(ggplot(as.data.frame(seur[[]][cells.to.keep& seur$sample==i, ]), aes(1, nCount_RNA)) + 
   geom_violin(fill="gray80") + theme_classic()+ theme(axis.title.x = element_blank())+
   geom_jitter(height = 0, width = 0.3, aes(col=nUMI.outlier.low)) +
   scale_color_manual(values=c("#00BFC4", "#F8766D"))+ggtitle(paste0("Total UMI counts per cell: ",i)))
}

for ( i in unique(seur$sample)){
  print(ggplot(as.data.frame(seur[[]][ cells.to.keep, ]), aes(1, nFeature_RNA)) + 
   geom_violin(fill="gray80") + theme_classic()+ theme(axis.title.x = element_blank())+
   geom_jitter(height = 0, width = 0.3, aes(col=nGene.outlier.low)) +
   scale_color_manual(values=c("#00BFC4", "#F8766D"))+ggtitle(paste0("Number of genes per cell: ",i)))
}

for ( i in unique(seur$sample)){
  print(ggplot(as.data.frame(seur[[]][cells.to.keep, ]), aes(1, perc.mito)) + 
   geom_violin(fill="gray80") + theme_classic()+ theme(axis.title.x = element_blank())+
   geom_jitter(height = 0, width = 0.3, aes(col=mito.outlier.high)) +
   scale_color_manual(values=c("#00BFC4", "#F8766D"))+ggtitle(paste0("% mito genes per cell: ",i)))
}

Filter the outlier cells

  seur_clean<-seur[, cells.to.keep]
  print(dim(seur))

[1] 28692 28284

  print(dim(seur_clean))

[1] 28692 26760

2.2) Get PCA outliers (Optional)

Get additional outliers, using multivariate outlier detection function of scater, based on PCA computed on QC metric data (uses the mvoutlier package)

sce=as.SingleCellExperiment(seur_clean)
##### Choose QC variables for authomatic outlier detection
selected_variables <- c("nCount_RNA", 
                        "nFeature_RNA", "perc.mito")
setdiff(selected_variables, colnames(colData(sce)))

character(0)

Detect outliers

outliers=c()
for ( i in unique(seur$sample)){
sce_sample=runColDataPCA(sce[,sce$sample==i], outliers=T, variables=selected_variables)
outliers=c(outliers,sce_sample$outlier)
}

sce$outliers=outliers
seur_clean <- AddMetaData(object = seur_clean, metadata = colData(sce)[colnames(seur_clean),"outliers"], col.name = "pca.outlier")
cat("Outliers:",sum(seur_clean$pca.outlier))

Outliers: 430

for ( i in unique(seur$sample)){
print(ggplot(as.data.frame(colData(sce)[sce$sample==i,]), aes(1, nCount_RNA)) + 
  geom_violin(fill="gray80") + theme_classic()+ theme(axis.title.x = element_blank())+
  geom_jitter(height = 0, width = 0.3, aes(col=outliers)) +
  scale_color_manual(values=c("#00BFC4", "#F8766D"))+ggtitle(paste0("Total UMI counts per cell: ",i)))
}

for ( i in unique(seur$sample)){
print(ggplot(as.data.frame(colData(sce)[sce$sample==i,]), aes(1, nFeature_RNA)) + 
  geom_violin(fill="gray80") + theme_classic()+ theme(axis.title.x = element_blank())+
  geom_jitter(height = 0, width = 0.3, aes(col=outliers)) +
  scale_color_manual(values=c("#00BFC4", "#F8766D"))+ggtitle(paste0("Number of genes per cell: ",i)))
}

  for ( i in unique(seur$sample)){
print(ggplot(as.data.frame(colData(sce)[sce$sample==i,]), aes(1, perc.mito)) + 
  geom_violin(fill="gray80") + theme_classic()+ theme(axis.title.x = element_blank())+
  geom_jitter(height = 0, width = 0.3, aes(col=outliers)) +
  scale_color_manual(values=c("#00BFC4", "#F8766D"))+ggtitle(paste0("% mito genes per cell: ",i)))
  }

As the number of additional PCA outliers is large, we keep those cells for now and will consider later if they need to be removed or not

### Remove the PCA outlier cells
# seur_clean<-seur_clean[, cells.to.keep]
# print(dim(seur_clean))

3) QC: GENES

To define a cutoff of lowly-abundant genes, we plot the distribution of log-means across all genes. The cutoff is placed in middle of the rectangular component of the graph before the peak.

thresholds<-c(0.005,0.005)
ave.counts=list()
for ( i in 1: length(unique(seur_clean$sample))) {
   ave.counts[[i]] <- rowMeans(as.matrix(GetAssayData(seur_clean, slot="counts")[,seur_clean$sample==unique(seur$sample)[i]]))
  hist(log10(ave.counts[[i]]), breaks=100, main=paste0("Histogram of mean UMI counts per gene: ",unique(seur$sample)[i]), col="grey80",
     xlab=expression(Log[10]~"mean count per gene"))
  abline(v=log10(thresholds[i]), col="blue", lwd=2, lty=2)
}

Number of genes to keep

[1] "KO1"

FALSE  TRUE 
16241 12451 
[1] "WT1"

FALSE  TRUE 
16471 12221 

Filter out the lowly-abundant genes that overlap between all samples

i=1
genes.filter=names(usegenes[[i]][! usegenes[[i]]])
for ( i in 2: length(unique(seur_clean$sample))) {
  genes.filter=intersect(genes.filter, names(usegenes[[i]][! usegenes[[i]]]))
}
usegenes.final=!rownames(seur_clean) %in% genes.filter
table(usegenes.final)

usegenes.final
FALSE  TRUE 
16021 12671 

seur_clean<-seur_clean[usegenes.final, ]

4) CALCULATE DOUBLET SCORE

We estimate doublet score per cell, using the scDblFinder package,which simulates artificial doublets from cell clusters. It is preferable to look for doublets separately for each sample

doublet.score<- scDblFinder::scDblFinder(as.SingleCellExperiment(seur_clean), samples="sample", BPPARAM=
BiocParallel::MulticoreParam(3),returnType="table")
seur_clean <- AddMetaData(object = seur_clean, metadata = doublet.score[colnames(seur_clean),"score"], col.name = "doublet.score")
seur_clean <- AddMetaData(object = seur_clean, metadata = doublet.score[colnames(seur_clean),"class"], col.name = "doublet.class")

table(seur_clean$doublet.class)

doublet singlet 
   3881   22879 

summary(seur_clean$doublet.score)

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
0.0000695 0.0078344 0.0356150 0.1956456 0.1480074 0.9999810 

5) NORMALIZATION, HVG DETECTION and PCA

##### Normalize data
seur_clean <- NormalizeData(seur_clean,verbose = F)
##### HVG detection 
seur_clean <- FindVariableFeatures(seur_clean,verbose=F)
##### Scale data per gene
seur_clean <- ScaleData(seur_clean,verbose=F)
##### PCA
seur_clean <- RunPCA(seur_clean, features =VariableFeatures(seur_clean) ,verbose=F)

Select PCs for downstream analysis of the dataset

Seurat provides a heuristic method to help us select PC components. It generates a ranking of PCs based on the percentage of variance explained by each of them

ElbowPlot(object = seur_clean,ndims =50)

Heatmap of the genes that drive each PC

DimHeatmap(seur_clean, dims = 1:30, cells = 5000, balanced = TRUE)

6) UMAP DIMENSIONALITY REDUCTION

Run Uniform Manifold Approximation and Projection (UMAP) dimensional reduction technique for visualisation of the data

### PC selection for downstream analysis
dims.use<-30

seur_clean <- RunUMAP(seur_clean, dims = 1:dims.use, verbose=F)

Visualise the UMAP plot, coloured by sample to check for batch effects

DimPlot(object = seur_clean, group.by = "sample",label=T, repel=T)

FeaturePlot(object = seur_clean, features = c("Ly6c2", "Ccr2"), split.by = "sample" )

The UMAP plot shows relatively good mixing by sample for most cells. Logically, monocytes are mainly present in Ccr2 WT cells, as visualized by Ccr2 and Ly6c2 expression.

7) BATCH CORRECTION

In case batch correction is necessary, it can be performed using the harmony package. Select a value for theta - diversity clustering penalty parameter (Default theta=2). Larger values of theta result in stronger integration, but can lead to over-correction)

theta.use<-1 # Here the theta parameter can be modified
seur_clean<-RunHarmony(seur_clean, group.by.vars="sample",theta =theta.use,  
                       plot_convergence = TRUE, reduction.save =paste0("harmonyTheta",theta.use),
                       reduction.key = paste0("harmonyTheta",theta.use,"_"), verbose=F) 

seur_clean <- RunUMAP(seur_clean,reduction =paste0("harmonyTheta",theta.use), 
                      dims = 1:dims.use,verbose =F,   reduction.name = paste0("umapHarmonyTheta",theta.use, "PC",dims.use), 
                      reduction.key = paste0("umapHarmonyTheta",theta.use, "PC",dims.use,"_"))  

### Visualise the harmony-corrected UMAP plot, coloured by sample to inspect if the batch effects were resolved
 DimPlot(object = seur_clean,reduction = paste0("umapHarmonyTheta",theta.use, "PC",dims.use), 
         group.by = "sample",label=T, repel=T)+ggtitle(paste0("Samples: UMAP Harmony theta= ",theta.use, " PC= ",dims.use))

As the the samples already show good mixing before batch correction, and the library preparation was done on the same day for both samples, we are going to continue without batch correction

Visualize QC metrics

FeaturePlot(object = seur_clean,
            features = c("nCount_RNA", "nFeature_RNA", "perc.mito","doublet.score"))

Visualize PCA outlier cells

DimPlot(object = seur_clean, group.by = "pca.outlier") 

Visualize QC metrics for pca outliers, and the remaining cells

FeaturePlot(seur_clean, split.by ="pca.outlier" ,features=c("nCount_RNA", "nFeature_RNA", "perc.mito", "Mki67","Stmn1"), 
keep.scale="feature", pt.size = 0.2)

The PCA outliers include cells with high percent mitochondrial genes, as well as cells with high UMI content and number of genes. Many of the latter express proliferation markers (Mki67, Stmn1), the active proliferation explaining their higher number of UMIs and genes. Therefore the PCA outliers will not be removed.

8) CLUSTERING

We run Louvain clustering with varying the clustering resolution between 0 and 2

seur_clean <- FindNeighbors(seur_clean, dims = 1:dims.use, 
                      graph.name = paste0("RNA_snn_PC",dims.use), verbose=F)
for ( i in seq(0,2, 0.25))
  seur_clean <- FindClusters(seur_clean, resolution = i, graph.name=paste0("RNA_snn_PC",dims.use), verbose=F)

If the batch correction is necessary for a specific dataset, clustering can be performed on the harmony-corrected PCA embeddings:

#seur_clean <- FindNeighbors(seur_clean, dims = 1:dims.use, reduction =paste0("harmonyTheta",theta.use), graph.name = paste0("RNA_snn_harmony_theta",theta.use, ".PC",dims.use), verbose=F)
#for ( i in seq(0,2, 0.25))
#  seur_clean <- FindClusters(seur_clean, resolution = i, graph.name=paste0("RNA_snn_harmony_theta",theta.use, ".PC",dims.use), verbose=F)

Choosing a clustering resolution value

Plot of a clustering tree showing the relationship between clustering at different resolution (using the clustree package).This plot allows us to see how are the clusters related to each other and which ones are stable across different resolutions

clustree(seur_clean, prefix = paste0("RNA_snn_PC",dims.use,"_res."))+
  ggtitle(paste("PC =",dims.use))

The `add` argument of `group_by()` is deprecated as of dplyr 1.0.0.
Please use the `.add` argument instead.
This warning is displayed once every 8 hours.
Call `lifecycle::last_warnings()` to see where this warning was generated.

The clusters are relatively stable until resolution of 1

Visualize the clusters for several resolution values on the UMAP plot

plot<-list()
for ( res in c(0.25,0.5, 0.75,1))
  plot[[as.character(res)]]<-DimPlot(seur_clean,label=T,repel=T, group.by = paste0("RNA_snn_PC",dims.use,"_res.",res))+
                                      ggtitle(paste("PC =",dims.use,"res=",res))
plot_grid(plotlist=plot)

9) ANNOTATING THE CLUSTERS

Let’s find differentially expressed (DE) genes per cluster for resolution=1.

res=1
Idents(seur_clean)=  paste0("RNA_snn_PC",dims.use,"_res.",res)
Idents(seur_clean)=  factor(Idents(seur_clean),levels = 0:(length(unique(Idents(seur_clean)))-1))
DEgenes<-FindAllMarkers(seur_clean,min.cells.group=2,pseudocount.use = 0.01, max.cells.per.ident = 1000,verbose = F)

Next,we visualise the expression of the top marker genes for each cluster in a heatmap

features.plot=DEgenes%>%filter(avg_logFC>0)%>% group_by(cluster)%>%top_n(n=-3,wt=p_val_adj)%>%pull(gene) ### top 3 markers
length(features.plot)

[1] 83

seur_clean=ScaleData(seur_clean, features = rownames(seur_clean)) ### Scale the data for all genes to be able to visualise their expression with DoHeatmap

DoHeatmap(seur_clean, features = features.plot, assay = "RNA", angle = 90, label =T, size=4) +
  scale_fill_gradient2(low = "blue", mid = "white",high = "red")+
  theme(legend.position = "bottom")%>% suppressMessages()

Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing scale.

Check the distribution of each cluster per sample

data<-as.data.frame(table(seur_clean$sample, Idents(seur_clean)))
colnames(data)=c("sample","cluster","Freq")
ggplot(data, aes(x=cluster, y=Freq, fill=sample))+geom_bar(stat = "identity")+ggtitle(paste0("Number of cells per cluster (resolution =",res,")"))+theme_classic()

Visualise the expression of marker genes on a UMAP plot

FeaturePlot(seur_clean,features=c("Ccr2","Ly6c2","C1qa","Sparc","P2ry12","Tgfbi","Ms4a7","Mki67","Gata2","Hdc","Cpa3","S100a9","Csf3r","Flt3","Xcr1","Cd209a","H2-DMb2","Ccr7","Ccr9","Ms4a1","Jchain", "Cd3e","Cd8a","Cd4","Foxp3","Klrb1c","Sox2","Plp1"))

seur_clean$nCount_RNA_log=log2(seur_clean$nCount_RNA)
seur_clean$nFeature_RNA_log=log2(seur_clean$nFeature_RNA)

Visualise the expression of marker genes using a DotPlot

DotPlot(seur_clean,features=c("Ptprc","Ccr2","Ly6c2","Plac8","C1qa","Sparc","Sall1","Tgfbi","Ms4a7","Mki67","Stmn1","Mcm3","Ifit3","Ly6g","Csf3r","S100a9","S100a8","Flt3","Xcr1","Cd209a","H2-DMb2","Napsa","Siglech","Ccr9","Ccr7","Cacnb3","Ccl22", "Cd3e","Cd8a","Cd4","Lef1","Tcf7","Foxp3","Il2ra","Tcrg-C1","Trdc","Klrb1c","Ncr1","Gnly", "Cd19","Ms4a1","Sdc1","Jchain","Ms4a2", "Gata2","Cd200r3","Cdh1","Hgf","Hdc","Mcpt1","Cpa3","Sox2","Gjc3","Plp1", "doublet.score","perc.mito","nCount_RNA_log","nFeature_RNA_log"))+RotatedAxis()

Assign cell types to the clusters, based on the marker gene expression

Idents(seur_clean) <- plyr::mapvalues(x = Idents(seur_clean), from = 0:25, 
     to =c("TAM 1","CD8 T cells","TAM 2","TAM 3","TAM 4","Artefact","NK cells","cDC2","Basophils","Treg cells","TAM Proliferating","migDC","Monocytes","cDC1","T cells Proliferating","CD4 T cells 1","pDC","CD4 T cells 2","TAM 5","Plasma cells","B cells","TAM/T doublets","Mast cells","Neutrophils","Oligodendrocytes","TAM/Basophil doublets"))

DimPlot(seur_clean,repel =T,label=T) 

DotPlot(seur_clean,features=c("Ptprc","Ccr2","Ly6c2","Plac8","C1qa","Sparc","Sall1","Tgfbi","Ms4a7","Mki67","Stmn1","Mcm3","Ifit3","Ly6g","Csf3r","S100a9","S100a8","Flt3","Xcr1","Cd209a","H2-DMb2","Napsa","Siglech","Ccr9","Ccr7","Cacnb3","Ccl22", "Cd3e","Cd8a","Cd4","Lef1","Tcf7","Foxp3","Il2ra","Tcrg-C1","Trdc","Klrb1c","Ncr1","Gnly", "Cd19","Ms4a1","Sdc1","Jchain","Ms4a2", "Gata2","Cd200r3","Cdh1","Hgf","Hdc","Mcpt1","Cpa3","Sox2","Gjc3","Plp1", "doublet.score","perc.mito","nCount_RNA_log","nFeature_RNA_log"))+RotatedAxis()

10) REMOVE DOUBLETS AND ARTEFACTS

Remove the doublets and artefacts - cluster with high mitochodrial content and the oligodendrocytes as the latter are probably dues to impurities during the CD45+ sorting

data<-as.data.frame(table(seur_clean$doublet.class, Idents(seur_clean)))
colnames(data)=c("doublet.class","cluster","Freq")
ggplot(data, aes(x=cluster, y=Freq, fill=doublet.class))+geom_bar(stat = "identity")+ggtitle(paste0("Number of cells per cluster (resolution =",res,")"))+theme_classic()+RotatedAxis()

cells.keep=WhichCells(seur_clean, idents =c("Artefact","TAM/T doublets","TAM/Basophil doublets","Oligodendrocytes"), invert = T )
cat(ncol(seur_clean)-length(cells.keep),"cells to be removed")

2350 cells to be removed

Remove cells and rerun HVG detection,PCA and UMAP

seur_artefacts_removed=seur_clean[,cells.keep]
seur_artefacts_removed <- FindVariableFeatures(seur_artefacts_removed,verbose=F)
seur_artefacts_removed <- ScaleData(seur_artefacts_removed,verbose=F)
seur_artefacts_removed <- RunPCA(seur_artefacts_removed, features =VariableFeatures(seur_artefacts_removed) ,verbose=F)

#### Select PCs for downstream analysis of the dataset
ElbowPlot(object = seur_artefacts_removed,ndims =50)

Heatmap of the genes that drive each PC

DimHeatmap(seur_artefacts_removed, dims = 1:30, cells = 5000, balanced = TRUE)

Run UMAP

dims.use<-30 ### Final PC selection for downstream analysis
seur_artefacts_removed <- RunUMAP(seur_artefacts_removed, dims = 1:dims.use, verbose=F)

DimPlot(object = seur_artefacts_removed, group.by = "sample",label=T, repel=T)+
  ggtitle(paste0("Samples: UMAP PC= ",dims.use))

Explore DE genes between cluster TAM5 and TAM1 through Volcano plot

contrast.name="TAM5vsTAM1"
DEG_TAM5vsTAM1=FindMarkers(seur_artefacts_removed,ident.1 = "TAM 5",ident.2 = "TAM 1",
                           pseudocount.use = 0.01, verbose = F)

DEG_TAM5vsTAM1$gene=rownames(DEG_TAM5vsTAM1)

ggplot(data=DEG_TAM5vsTAM1, aes(x= avg_logFC, y = -log10(p_val_adj))) +
    geom_point(alpha = 1, size=2, 
               aes(col = p_val_adj<1e-20 & (avg_logFC>1 | avg_logFC< -1))) +
    scale_color_manual(values = c("dark grey", "brown1")) +
    geom_text_repel(data=subset(DEG_TAM5vsTAM1,  
                    p_val_adj<1e-100 & (avg_logFC>1 | avg_logFC< -1)),aes(label = gene), size =4)+
    theme_classic()+theme(legend.position = "none")+ggtitle(contrast.name)

Gene Ontology enrichment with Metascape of the genes upregulated in cluster TAM5 compared to cluster TAM1

Select top upregulated DEG (as those genes are specific for the cluster of interest). It is best to have at least 20-50 genes for GO enrichment.

p_tresh=1e-20
logFC_thresh=1
selected.genes=DEG_TAM5vsTAM1[DEG_TAM5vsTAM1$p_val_adj<p_tresh & DEG_TAM5vsTAM1$avg_logFC>logFC_thresh,"gene"]
length(selected.genes)

[1] 104

Export gene names as a csv file for metascape

write.table(selected.genes, file=paste0("Upregulated_genes_",contrast.name,"_p_val_adj.",p_tresh,"_logFC.",logFC_thresh,".csv"), row.names = F, col.names = F, sep=",")

Run GO enrichment in Metascape (https://metascape.org/):

Input the .csv file and do Express analysis with “Input as”= M.musculus and “analysis as”= M.musculus => Express Analysis. Save the “Gene List Report Excel Sheets”

Visualization of the enriched pathways on the Volcano plot:

GO.metascape=xlsx::read.xlsx( paste0("Metascape_Upregulated_genes_",contrast.name,"_p_val_adj.",p_tresh,"_logFC.",logFC_thresh,".xlsx"),2 )
#Select only the main (summary) pathways
GO.metascape$group=sapply(strsplit(GO.metascape$GroupID, "_"), "[[",2)

Top GO enriched pathway:

index=1
gset.symbols=GO.metascape[index,"Symbols"]
gset.symbols=strsplit(gset.symbols, ",")[[1]]
GOterm=GO.metascape[index,"Description"]
ggplot(data=DEG_TAM5vsTAM1, aes(x= avg_logFC, y = -log10(p_val_adj))) +
    geom_point(alpha = 1, size=2,
               aes(color=gene %in% gset.symbols)) +
    scale_color_manual(values = c("dark grey", "brown1")) +
    geom_text_repel(data=subset(DEG_TAM5vsTAM1,  
                   gene %in% gset.symbols),aes(label = gene), size =4)+
    theme_classic()+theme(legend.position = "none")+ggtitle(paste0( contrast.name, ": ", GOterm, " pathway"))

“response to hypoxia” pathway:

GOterm="response to hypoxia"
gset.symbols=GO.metascape[GO.metascape$Description==GOterm,"Symbols"]
gset.symbols=strsplit(gset.symbols, ",")[[1]]
ggplot(data=DEG_TAM5vsTAM1, aes(x= avg_logFC, y = -log10(p_val_adj))) +
    geom_point(alpha = 1, size=2,
               aes(color=gene %in% gset.symbols)) +
    scale_color_manual(values = c("dark grey", "brown1")) +
    geom_text_repel(data=subset(DEG_TAM5vsTAM1,  
                   gene %in% gset.symbols),aes(label = gene), size =4)+
    theme_classic()+theme(legend.position = "none")+ggtitle(paste0( contrast.name, ": ", GOterm, " pathway"))

DimPlot(object = seur_artefacts_removed,label=T, repel=T)+
  ggtitle(paste0("Cell annotation: UMAP PC= ",dims.use))

Save object

saveRDS(seur_clean, file="Citeseq_mouse_GBM.seuratObj.rds")

sessionInfo()

R version 4.0.3 (2020-10-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.4 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0

locale:
 [1] LC_CTYPE=en_US.UTF-8          LC_NUMERIC=C                  LC_TIME=en_US.UTF-8           LC_COLLATE=en_US.UTF-8       
 [5] LC_MONETARY=en_US.UTF-8       LC_MESSAGES=en_US.UTF-8       LC_PAPER=en_US.UTF-8          LC_NAME=en_US.UTF-8          
 [9] LC_ADDRESS=en_US.UTF-8        LC_TELEPHONE=en_US.UTF-8      LC_MEASUREMENT=en_US.UTF-8    LC_IDENTIFICATION=en_US.UTF-8

attached base packages:
 [1] grid      parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] ggrepel_0.9.1               scDblFinder_1.4.0           cowplot_1.1.1               clustree_0.4.3             
 [5] ggraph_2.0.4                harmony_1.0                 Rcpp_1.0.6                  mvoutlier_2.0.9            
 [9] sgeostat_1.0-27             VennDiagram_1.6.20          futile.logger_1.4.3         dplyr_1.0.3                
[13] scater_1.18.3               ggplot2_3.3.4               SingleCellExperiment_1.12.0 SummarizedExperiment_1.20.0
[17] Biobase_2.50.0              GenomicRanges_1.42.0        GenomeInfoDb_1.26.7         IRanges_2.24.1             
[21] S4Vectors_0.28.1            BiocGenerics_0.36.1         MatrixGenerics_1.2.1        matrixStats_0.60.0         
[25] Seurat_3.2.3               

loaded via a namespace (and not attached):
  [1] scattermore_0.7           prabclus_2.3-2            GGally_2.1.2              tidyr_1.1.2               knitr_1.30               
  [6] irlba_2.3.3               DelayedArray_0.16.3       data.table_1.13.6         rpart_4.1-15              RCurl_1.98-1.4           
 [11] generics_0.1.0            lambda.r_1.2.4            RANN_2.6.1                proxy_0.4-24              future_1.21.0            
 [16] spatstat.data_2.1-0       httpuv_1.5.5              assertthat_0.2.1          viridis_0.5.1             xfun_0.20                
 [21] rJava_0.9-13              hms_1.0.0                 evaluate_0.14             promises_1.1.1            fansi_0.4.2              
 [26] DEoptimR_1.0-8            readxl_1.3.1              igraph_1.2.6              DBI_1.1.1                 htmlwidgets_1.5.3        
 [31] reshape_0.8.8             purrr_0.3.4               ellipsis_0.3.1            RSpectra_0.16-0           ks_1.13.1                
 [36] backports_1.2.1           sROC_0.1-2                deldir_1.0-6              sparseMatrixStats_1.2.0   vctrs_0.3.6              
 [41] ROCR_1.0-11               abind_1.4-5               withr_2.4.0               ggforce_0.3.2             robustbase_0.93-7        
 [46] checkmate_2.0.0           vcd_1.4-8                 sctransform_0.3.2         scran_1.18.3              mclust_5.4.7             
 [51] goftest_1.2-2             cluster_2.1.0             lazyeval_0.2.2            laeken_0.5.1              crayon_1.4.0             
 [56] labeling_0.4.2            edgeR_3.32.1              pkgconfig_2.0.3           zCompositions_1.3.4       tweenr_1.0.1             
 [61] nlme_3.1-149              vipor_0.4.5               nnet_7.3-14               rlang_0.4.10              globals_0.14.0           
 [66] diptest_0.75-7            lifecycle_0.2.0           pls_2.7-3                 miniUI_0.1.1.1            rsvd_1.0.3               
 [71] cellranger_1.1.0          polyclip_1.10-0           lmtest_0.9-38             Matrix_1.3-2              carData_3.0-4            
 [76] boot_1.3-25               zoo_1.8-8                 beeswarm_0.2.3            ggridges_0.5.3            png_0.1-7                
 [81] viridisLite_0.3.0         bitops_1.0-7              KernSmooth_2.23-17        DelayedMatrixStats_1.12.2 stringr_1.4.0            
 [86] parallelly_1.23.0         beachmat_2.6.4            scales_1.1.1              magrittr_2.0.1            plyr_1.8.6               
 [91] ica_1.0-2                 zlibbioc_1.36.0           compiler_4.0.3            hdrcde_3.4                dqrng_0.2.1              
 [96] RColorBrewer_1.1-2        rrcov_1.5-5               fitdistrplus_1.1-3        cli_2.2.0                 XVector_0.30.0           
[101] listenv_0.8.0             patchwork_1.1.1           pbapply_1.4-3             formatR_1.7               MASS_7.3-53              
[106] mgcv_1.8-33               tidyselect_1.1.0          stringi_1.5.3             forcats_0.5.0             yaml_2.2.1               
[111] BiocSingular_1.6.0        locfit_1.5-9.4            tools_4.0.3               future.apply_1.7.0        rio_0.5.16               
[116] rstudioapi_0.13           bluster_1.0.0             foreign_0.8-79            gridExtra_2.3             farver_2.0.3             
[121] Rtsne_0.15                digest_0.6.27             shiny_1.5.0               pracma_2.3.3              fpc_2.2-9                
[126] car_3.0-10                xlsx_0.6.5                scuttle_1.0.4             later_1.1.0.1             RcppAnnoy_0.0.18         
[131] fda_5.1.9                 httr_1.4.2                kernlab_0.9-29            colorspace_2.0-0          tensor_1.5               
[136] rainbow_3.6               ranger_0.12.1             reticulate_1.18           truncnorm_1.0-8           splines_4.0.3            
[141] uwot_0.1.10               statmod_1.4.35            spatstat.utils_2.2-0      graphlayouts_0.7.1        xlsxjars_0.6.1           
[146] sp_1.4-5                  xgboost_1.3.2.1           flexmix_2.3-17            plotly_4.9.3              xtable_1.8-4             
[151] jsonlite_1.7.2            fds_1.8                   futile.options_1.0.1      spatstat_1.64-1           tidygraph_1.2.0          
[156] modeltools_0.2-23         R6_2.5.0                  NADA_1.6-1.1              pillar_1.4.7              htmltools_0.5.2          
[161] mime_0.9                  glue_1.4.2                fastmap_1.1.0             VIM_6.1.0                 BiocParallel_1.24.1      
[166] BiocNeighbors_1.8.2       cvTools_0.3.2             class_7.3-17              codetools_0.2-16          pcaPP_1.9-74             
[171] mvtnorm_1.1-1             lattice_0.20-41           tibble_3.0.5              curl_4.3                  ggbeeswarm_0.6.0         
[176] leiden_0.3.6              zip_2.1.1                 openxlsx_4.2.3            survival_3.2-7            limma_3.46.0             
[181] rmarkdown_2.6             munsell_0.5.0             e1071_1.7-6               GenomeInfoDbData_1.2.4    haven_2.3.1              
[186] reshape2_1.4.4            gtable_0.3.0              robCompositions_2.3.0